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I review the present status of nonperturbative studies on the nature of the fi- 
nite temperature QCD transition on the lattice, with the special attention to the 
determination of the order of the transition. 

1 Introduction 

Unlike in the case of the usual formulation of renormalizable field theories in 
continuum space-time, where the UV cutoff is introduced through the pertur- 
bation theory, the lattice formulation 1 uses a finite lattice spacing, a, which 
acts as a UV cut-off independent of the validity of the perturbative expansions. 
This enables us to study nonperturbative properties of the theory, which are 
essential in the finite temperature physics of the quantum chromodynamics 
(QCD). 

The continuum theory is defined in the limit of vanishing lattice spacing, 
a = 0, and infinite lattice volume. Before taking the continuum limit, the 
theory is finite and well defined, so that we can perform numerical simulations 
on super computers. Extrapolating the results to the continuum limit, we 
are discussing, for example, the hadron spectrum and decay constants to the 
accuracy of about 10%. 

In this paper, I review the present status of numerical lattice calculations 
for the finite temperature QCD transition, with the special attention to the 
determination of the order of the transition taking the effects of dynamical 
quarks into account. For recent reviews on finite temperature QCD on the 
lattice, see Refs. 2 and 3. 
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In nature, we know six flavors of quarks: u, d, s, c, b, and t. The lightest 
two quarks, u and d, are much lighter than the relevant energy scales whose 
typical magnitude can be given by the critical temperature T c : m u , <C T c 
(in the units h = c = kg = 1). On the other hand, the last three quarks, 
c, b, and t, are sufficiently heavy and are expected to play no appreciable 
roles to the thermal processes near T c . In this review, I concentrate on the 
physically interesting cases of two and three flavors of degenerate light quarks 
(Np = 2 and 3) and the more realistic case of two light and one heavy quarks 
(Np = 2 + 1). The case Np = 2 corresponds to the case that the third quark, 
the s quark, is much heavier than the relevant energy scale: m s ^> T c . The case 
Np = 3 corresponds to the case m s <C T c . In the following, we will see that 
the nature of the transition in the chiral limit (the limit of vanishing quark 
masses) for Np = 2 is different from that for Np = 3. Because m s ~ 150 
- 200 MeV is just of the same order of magnitude as the expected values of 
T c , this difference means that, in order to make a reliable prediction for the 
real world, we have to fine-tune the value of m s in the more realistic case 
Np = 2 + 1. Unfortunately, the lattice simulations performed so far do not 
have that accuracy. 

In the next section, I give a short introduction to the lattice formulation 
of QCD. In Sect. 3, I review recent studies on the QCD transition in the chiral 
limit. I then discuss the influence of the s quark in Sect. 4. A brief summary 
is given in Sect. 5. 

2 QCD on the lattice 

On a 4-dimensional Euclidian lattice, the gluon fields of QCD are provided by 
3x3 complex matrices, U Xtfl £ SU(3), living on the "links" connecting two 
neighboring sites, x and x + /}, where ft is the lattice unit vector in the fi 
direction. Near the continuum limit, a ~ 0, the "link variable" U Xtfl is related 
to the conventional gluon field A^(x) by 

U x>ll = exp[igaA^(x)}, (1) 

where g is the gauge coupling constant. As in the continuum limit, quark fields 
on the lattice are given by Grassmann variables living on the sites. 

The action on the lattice is chosen such that the continuum action is 
recovered when we let a — ► in the action. We also want to keep as much 
symmetries of the continuum action as possible on the lattice. 

The minimal choice for the gluon action is the standard one-plaquette 
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action by Wilson, 4 

Sgluon = (3 ^2 ^R-eTWp/a, (2) 

plaquettes 

where U p i aq is the ordered product of four JX^'s winding around a minimal 
square, the "plaquette", and (3 = 6/g 2 . Due to the quantum fluctuations, the 
lattice spacing a is a decreasing function of (3 described by the renormalization 
group beta function. Therefore, the continuum limit corresponds to the limit 

/3 = oo. 

On the lattice, however, we have a freedom to introduce less local terms 
without affecting the continuum limit. By adjusting these additional coupling 
parameters, it is possible to accelerate the approach to the continuum limit. 
Such actions are called "improved actions" . Importance of improvement is 
much stressed recently in many studies of lattice gauge theories. I will discuss 
an example for finite temperature physics later. 

For the quark part, it is known that naive lattice discretizations of the 
continuum action have the "doubler problem" : they necessarily contain addi- 
tional fermionic degrees of freedom which survive in the continuum limit. It is 
inevitable to adopt a slightly complicated formulation on the lattice, abandon- 
ing some of simple features of the continuum theory. 5 Two conventional choices 
are the staggered and the Wilson fermions. 

In the formulation of staggered quarks, 6 we introduce one component 
Grassmann field \ on each site. The 4 component Dirac field "3/ is constructed 
collecting x fields distributed on hypercube. Because a hypercube has 2 4 sites, 
we end up with 4 flavors of Dirac fields. These 4 flavors are degenerate also 
in the continuum limit. Due to a flavor mixing interaction among these 4 fla- 
vors, the flavor/chiral symmetry SV(Np)i x SV(Np)n x U(l) of massless QCD 
in the continuum limit b is explicitly broken down to \J(Np/4)i x ~U(Np /4)r 
at a > 0. But, because at least a part of the chiral symmetry is preserved, 
the location of the chiral limit is protected by this symmetry in the coupling 
parameter space of the lattice action. This makes the choice of simulation pa- 
rameters easier. A big flaw of this formulation is that Np must be a multiple 
of 4. For the most interesting cases Np = 2 and 3, a usual trick is to modify 
the power of the fermionic determinant in the numerical path-integration by 
hand. This necessarily makes the action non-local, that sometimes poses con- 
ceptually difficult problems, as I discuss later. Practically, numerical results 
of hadron spectrum seem to suggest that the major lattice artifacts become 
quickly small at large (3. 

6 The classical symmetry, V(Np )l x V(Np )r, is broken down to this symmetry due to 
the quantum 11^(1) anomaly. 
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Another conventional choice for the quark action is the Wilson fermion 
action. 7 In this formulation, we introduce four component Dirac field "J" on 
each site. Therefore the flavor symmetry is manifest also on the lattice. In 
order to avoid the doubler problem, however, we have to introduce the "Wilson 
term" to the action, which vanishes when we let a — ► naively in the action but 
explicitly breaks the chiral symmetry at a > 0. The numerical studies show 
that the effects of this chiral violation are small and are well under control. 

In the continuum limit, we expect that both Wilson and staggered quarks 
give the same results, at least for Np = 4x integer. Off the continuum limit, 
they are in general different, and the difference provides us with a measure of 
the quality of predictions from the lattice. 

3 Finite temperature chiral transition in QCD 

At low temperatures, QCD has three characteristic properties: confinement of 
quarks and gluons, spontaneous breakdown of the chiral symmetry, and the 
asymptotic freedom. The first two are IR properties and may be modified by 
temperature. Actually, the study of lattice QCD has shown that there exists 
a finite temperature transition between the low temperature hadronic phase 
where quarks are confined and the chiral symmetry is spontaneously broken, 
and the high temperature quark gluon plasma (QGP) phase where quarks are 
deconfined and the chiral symmetry is restored. Lattice simulations further 
suggest that the confinement-deconfinement transition and the recovery of the 
chiral transition take place at the same temperature. 

In QCD we have several adjustable parameters such as quark masses and 
the number of flavors. Phenomenological interests are, of course, concentrated 
on the nature of the transition in a particular case where the quark masses have 
their physical values. This, however, requires a rather delicate fine-tuning of 
parameters in a lattice simulation. Therefore, it is important to make clear 
the nature of the transition in limiting cases. 

One of these interesting limiting cases of QCD is the limit of the SU(3) pure 
gauge theory, that corresponds to the case of infinitely heavy quarks, m q ^> T c . 
In this limit, QCD has a global Z(3) symmetry, which is spontaneously broken 
when quarks are deconfined. This suggests that the deconfining transition of 
the SU(3) gauge theory belongs to the same universality class as the order- 
disorder transition of three dimensional Z(3) spin models such as the 3-state 
Potts model. Because the transition is of first order in these spin models, we 
expect that the deconfinement transition is also of first order. Precise lattice 
studies applying the finite size scaling analysis confirmed this expectation. 8 ' 9 

The success of the universality argument in the pure gauge limit encour- 
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Table 1: Critical exponents of the three dimensional 0(4) 14 and 0(2) 15 Heisenberg models, 
and Np — 2 QCD with staggered quarks. 16 





0(4) 


0(2) 


Karsch-Laermann 


1/(36 


0.537(7) 


0.602(2) 


0.77(14) (3 C 


1/6 


0.2061(9) 


0.2072(3) 


0.21 - 0.26 chiral cumulant 


1 - 1/6 


0.7939(9) 


0.7928(3) 


0.79(4) X m 


(1-/3)/ /36 


0.331(7) 


0.394(2) 


0.65(7) X t 


a/136 


-0.13(3) 


-0.003(4) 


-0.07 - +0.34 specific heat 



ages us to apply a similar idea to more realistic cases with dynamical quarks. 
Here, the Z(3) symmetry of the pure gauge theory is explicitly broken by the 
quark mass term. Therefore, we cannot use this symmetry for universality 
argument. In the massless limit (chiral limit), however, we have the sponta- 
neously broken chiral symmetry that is expected to recover at T > T c . Based 
on a study of an effective a model which respects the chiral symmetry of QCD, 
Pisarski and Wilczek 10 discussed that the transition in the chiral limit (chiral 
transition) of QCD is of first order for Np > 3. This was confirmed by numer- 
ical simulations both with staggered 11 and with Wilson quarks. 12 (See Ref. 2 
for a more complete list of references.) 

For two flavors, theoretical studies do not have a definite prediction about 
the order of the chiral transition: the transition is either of second order or 
of first order depending on the strength of the LU(1) anomaly term at the 
transition temperature. 10 In case that the chiral transition is of second order, 
the universality argument suggests that it will belong to the same universality 
class as the order-disorder transition of the three dimensional 0(4) Heisenberg 
model, 13 which shows well-studied critical behaviors. This provides us with 
several useful scaling properties among physical quantities around the chiral 
transition that allow us to make clear the order of the transition by numerical 
simulations. 

Scaling properties are described in terms of critical exponents. Consider a 
spin model at temperature T near the transition temperature T c and at small 
external magnetic field h. The exponents a, (3, 7, and 6 are defined for the 
specific heat C , the magnetization M , and the magnetic susceptibility \ by 



C(t,h = 0) ~ \t\- a 




(3) 


M(t,h = 0) ~ \tf 


(t<0) 


(4) 


x (t,h = 0) ~ t' 1 


(t>0) 


(5) 


M(t = 0,h) ~ h 1/s 




(6) 



where t = [T — T c (h = Q)]/T c (h = Q) is the reduced temperature. These critical 
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exponents satisfy the following two scaling relations: a = 2 — (3(6 + 1) and 
7 = (3(6 — 1), so that only two of the four exponents are independent. We also 
expect that magnetization M near the transition point can be described by a 
single scaling function: 

M = h 1 / s f(t/h 1 ^ s ). (7) 

Some 0(4) exponents are listed in Table 1. 

In QCD, h corresponds to the quark mass and M corresponds to the chiral 
condensate (the chiral order parameter) ("I"]/}. On a lattice with the linear 
lattice size N t in the time direction, temperature is given by T = l/[N t a((3)], 
where the lattice spacing a is a decreasing function of (3, as discussed before. In 
simulations with N t fied, larger (3 corresponds to larger T. The continuum limit 
at a given value of T is achieved by letting N t — ► oo and a — ► 0, simultaneously. 

The conjecture of the 0(4) scaling has been tested first for staggered quarks 
by F. Karsch and E. Laermann on an N t = 4 lattice. 16 The last two columns 
in Table 1 show their numerical results for critical exponents. The result 
for 6 is consistent with the 0(4) value, while the results for a and (3 are in 
disagreement with the 0(4) exponents. It was argued that the discrepancies 
of a and (3 might be caused by the small lattice size simulated. C. DeTar 
tested the scaling relation (7) using data for (tyty) from two flavor staggered 
quarks 17 and found that data are roughly consistent with the 0(4) saling 
(and also with the 0(2) scaling, discussed below). With staggered quarks, 
however, several caveats are in order because Np = 2 staggered quarks are 
realized on the lattice by introducing a fractional exponent to the fermionic 
determinant of the Np = 4 staggered quarks. Therefore, (i) the symmetry 
in the chiral limit on a lattice with finite lattice spacing is not the 0(4) but 
0(2) of Np = 4 staggered quarks, and (ii) the action is not local. The correct 
continuum chiral limit with the 0(4) symmetry will be obtained only when we 
first take the continuum limit (3 — ► oo and then take the chiral limit. Choosing 
a too small m q compared with the lattice spacing may lead us either to wrong 
0(2) exponents, or to some non-universal behavior due to the lack of locality. 
In order to conclude the 0(4) scaling among these possibilities, we clearly 
need data with better accuracy on larger lattices. At least, the appearance of 
non-trivial exponents strongly suggests that the transition is of second order, 
although the lattice artifact may be large as we cannot see the clear 0(4) 
scaling. 

In the case of Wilson quarks, the lattice action lacks chiral symmetry 
due to the Wilson term. This introduces possible 0(a) effects to physical 
quantities and relations characteristic for chiral symmetry. However, when 
we are close enough to the continuum limit, we expect that the symmetry 
breaking effects become sufficiently weak so that we see the 0(4) scaling when 
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Figure 1: (a) Chiral condensate as a function of 2m q a for Wilson quarks with a RG improved 
action on an 8 3 X 4 lattice. 18 (b) Best fit for the scaling function with 0(4) exponents. The 
plot contains all data of Fig. 1(a) within the range < 2m q a < 0.9 and f} < 2.0. 



the chiral transition is of second order. 2 With the standard action, we found 
that the effects of the Wilson term are not negligible on lattices available with 
the present power of computers. 3 We therefore applied an improved action to 
study this issue. 18 We first confirmed that the lattice artifact from the Wilson 
term is small with our improved action. We then studied the scaling function 
in (7) with the identification h = 2m q a and t = (3 — f3 ct , where (3 ct is the 
chiral transition point. Fixing the exponents to 0(4) values given in Table 1, 
we adjust (3 ct to obtain the best fit. Our result is shown in Fig. 1. We find 
that the scaling ansatz works remarkably well. This result is quite nontrivial 
because a slight shift of the exponents to other values suggested, for example, 
from the meanfield theory makes the fit apparently worse. 



4 Influence of the strange quark 

In the previous section, we have seen that the finite temperature transition 
is of second order in the limit of massless quarks in Np = 2 QCD, while it 
is of first order for Np > 3. Off the chiral limit, the first order transition 
smoothens into a crossover at sufficiently large m q . Therefore, the nature of 
the transition sensitively depends on Np and m q . This means that, in order 
to study the nature of the transition in the real world, we should include the 
s quark properly whose mass m s is of the same order of magnitude as the 
transition temperature T c ~ 100 - 200 MeV. 

Following Brown et al.} 1 we summarize in Fig. 2 what we expect about 



7 



m. 





N F =2 


pure ga 


to 




\ 1st 


nd 




V>rder 


o 
Cu 


crossover 




CD 













1st 
order 











m 



ud 



Figure 2: Map of expected nature of the QCD transition for Np = 2 + 1 QCD as a function 
of the u and d quark mass m u d and the s quark mass m s . 



the nature of the finite temperature transition as a function of quark masses, 
neglecting the mass difference among u and d quarks {Np = 2 + 1). When all 
quarks are heavy, the transition is of first order as observed in the SU(3) pure 
gauge theory. The limit m s = oo corresponds to the case Np = 2 discussed in 
Sect. 3 where we found strong evidence for second order transition at m u d = 0. 
When the chiral transition is of second order, we expect that it turns into an 
analytic crossover at nonvanishing quark masses. For m u d = m s (Np = 3), the 
transition is of first order in the chiral limit. Therefore, on the axis m u d = 0, 
we have a tricritical point to* where the second order transition at large m s 
turns into first order 13 . For m s > to*, the second order transition line follows 
the m u d = axis and, for m s < m*, is suggested 19 to deviate from the vertical 
axis according to m u d oc (to* — to s ) 5 / 2 . 

Our main goal is to determine the position of the physical point in this 
map. 
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Np= 2 pure gauge 




Figure 3: Same as Fig. 2 with the results of simulations with Wilson quarks using the 
standard action. 20 First order signals are observed at the points marked by filled circles, 
while no clear two state signals are found at the points represented by the open circles. The 
values of quark mass in physical units are computed using a -1 ~ 0.8 GeV for /3 < 4.7 and 
a' 1 ~ 1.0(1.8) GeV for /3 = 5.0(5.5) determined by m p (T = 0) = 770 MeV. The real world 
corresponds to the point marked by the star. 



Using staggered quarks, the Columbia group studied this issue five years 
ago on JVj = 4 lattices. 11 For the degenerate Np = 3 case, m u d = m s = m q , 
they found a first order signal for m q a = 0.025 at (3 = 5.132. For Np = 2 + 1, 
they obtained a time history suggesting a crossover for m u da = 0.025, m s a = 
0.1 at (3 = 5.171. Their study of hadron spectrum at this point gives the value 
of nix/nip smaller than the experimental value, suggesting that their m s is 
smaller than its physical value. At the same time, their large m w /m p suggests 
that their m u d is larger than the physical value. This implies that the physical 
point is located in the crossover region unless the second order transition line, 
which has a sharp m u d dependence near to* (cf. Fig. 2), crosses between the 
physical point and the simulation point. 
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Recently, we studied this issue with Wilson quarks using the standard 
action on N t = 4 lattices. 20 For N F = 3, we performed simulations increasing 
the quark mass at the transition temperature and found that the first order 
signal observed in the chiral limit 12 persists for the cases corresponding to 
m q < 140 MeV, while no clear two state signals are observed for m q > 250 
MeV, where the physical s quark mass, giving = 1.02 GeV, is about 150 
MeV with our normalization of m q . For Np = 2 + 1, we observed first order 
signals for m u d ~ at both m s ~ 150 and 400 MeV. Our study of hadron 
spectroscopy for Np = 2 + 1 shows that ~ 1.03(5) GeV at the simulation 
point with m s ~ 150 MeV, verifying that this simulation point is very close to 
the physical point in this sense. Our results on the nature of the transition are 
summarized in Fig. 3. The physical point clearly falls in the first order region. 

Although both staggered and Wilson simulations give a phase structure 
qualitatively consistent with Fig. 2, Wilson quarks tend to give larger values for 
critical quark masses (measured by m,^,/m p etc.) than those with staggered 
quarks. This leads to the difference in the conclusions about the location 
of the physical point, in Fig. 2. On the other hand, both of these studies 
discuss that the deviation from the continuum limit is large at N t = 4 where 
these simulations are done: sizable deviations from the experimental values are 
observed in several physical quantities. We should certainly make a calculation 
at larger N t or with an improved action in order to draw a definite conclusion 
about the nature of the QCD transition in the real world. At present, increasing 
N t is quite painful. Therefore, we began to study these issues applying the 
improved action which was successful for Np = 2. Some preliminary results 
were reported, 18 so far with no conclusion about the values for critical quark 
masses. 

5 Conclusions 

Results from lattice QCD mean that they are derived from the first principles 
of the theory of strong interaction. They, however, require careful extrapo- 
lations to the continuum limit. Numerical simulations of lattice QCD show 
strong evidences that the QCD transition in the massless quark limit is of 
first order for Np > 3, while, for Np = 2, it is of second order universal to 
the order-disorder transition of the three dimensional 0(4) Heisenberg model. 
These results are in accordance with the theoretical expectations based on the 
universality argument. In order to make a prediction to the real world, we 
have to fine-tune the s quark mass to its physical value. Unfortunately, the 
present lattice simulations do not have that accuracy yet. We hope to solve 
this problem by new simulations applying improved actions, so that we may 



10 



soon be able to make more definite predictions about the natute of the finite 
temperature QCD transition at the early Universe. 
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